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A shadow is an exact solution to a chaotic system of equations that remains close to a numerically 
computed solution for a long time, ending in a glitch. We study the distribution of shadow durations 
at low dimension and how shadow durations scale as dimension increases up to 300 in a slightly 
simplified gravitational n-body system. We find that "softened" systems are shadowable for many 
tens of crossing times even for large n, while in an "unsoftened" system each particle encounters 
glitches independently as a Poisson process, giving shadow durations that scale as 1/n. 
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The behaviour of dynamical systems is often studied 
using numerical techniques. A source of error in such 
studies arises because dynamical systems often display 
sensitive dependence on initial conditions: two solutions 
whose initial conditions differ by an arbitrarily small 
amount generally diverge exponentially away from each 
other. Since numerical methods introduce errors, it is 
virtually guaranteed that a numerically computed solu- 
tion diverges exponentially away from the exact solution 
with the same initial conditions. This remains true even 
if integrals of motion such as energy and momentum are 
conserved to arbitrary precision. Although this fact is 
widely known, its impact is not well understood. 

Fortunately, most studies of dynamical systems do not 
aim to predict the precise evolution of a particular choice 
of initial conditions. Instead, the dynamics of the system 
is sampled in order to study its general behaviour. In 
such cases, we typically choose initial conditions from a 
random distribution and would be happy if our numerical 
solution exhibited behaviour typical of any valid choice of 
initial conditions from our distribution. In particular, we 
may be satisfied if our numerical solution closely follows 
some exact solution whose initial conditions are close to 
those that we chose. 

The study of shadowing provides just such a property: 
a shadow is an exact solution to a given set of equations 
that remains close to a numerically computed solution 
of the same set of equations for a nontrivial duration of 
time. Although not all numerical simulations are likely 
to be shadowable ||, |^, the existence of a shadow is 
a strong property: it asserts that a numerical solution 
can be viewed as an experimental observation of an ex- 
act solution. As such, within the "observational" error, 
the dynamics observed in a numerical solution that has 
a shadow represent the dynamics of an exact solution. 
There are only two remaining questions (both beyond 
the scope of this paper): (1) Whether the mathematical 
model being simulated accurately reflects the system be- 
ing studied. This is certainly not the case for systems 
such as the weather (and thus shadowing is probably an 
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inappropriate measure of error), but can be assumed to 
be the case for others, such as the unsoftened gravita- 
tional n-body problem. (2) Whether shadows are typical 
of exact solutions chosen at random. Simple examples 
exist of shadows that are atypical ||, ||, ||], although it 
seems unlikely that atypical shadows are common — lest 
the numerical solutions we compute would be commonly 
atypical as well p9[ . 

Shadows of "noisy" trajectories were first proved to ex- 
ist in hyperbolic systems [Q, ||. Such systems possess a 
phase space which is spanned locally by an independent 
set of stable and unstable directions that remain con- 
sistent under the evolution of the system. Few chaotic 
dynamical systems are globally hyperbolic, but for those 
that are nearly hyperbolic over finite time intervals, the 
existence of finite-duration shadows can often be inferred 
lH, |l^ or rigorously proved ^, |l2|. If a shadow is 
viewed as a measure of error of a numerical solution, 
then the relevent measures are the phase-space distance 
between corresponding points on the "noisy" and exact 
trajectories (smaller is better), and the duration over 
which they remain close together (longer is better). Gen- 
erally, the smaller the local error in the trajectory, and 
the more hyperbolic the trajectory, the closer and longer 
the shadow |]lOj . In |[ |j , scaling laws were developed 
for the expected duration of shadows, based upon the 
distance from zero of finite-time Lyapunov exponents of 
the system. If an exponent fluctuates about zero over 
a trajectory, it represents a tangential direction that is 
uncertain between stable and unstable. This effect, also 
called unstable dimension variability , is the cause of 
unshadowability in chaotic systems. If a numerical tra- 
jectory is unshadowable, then it is possible (though far 
from certain) that statistical quantities associated with 
the numerical trajectory can have significant bias ^, |^, ^. 

To date, most studies of shadowability of chaotic sys- 
tems have focussed on maps or ODEs with only a few 
dimensions, and concern has been expressed that fiuc- 
tuating Lyapunov exponents may be common in high- 
dimensional systems H, ^. Furthermore, since Hamil- 
tonian systems are not globally hyperbolic, there is no 
reason to expect them to be easily shadowable. In 
this paper, we demonstrate that trajectories of at least 
one commonly integrated continuous Hamiltonian system 



2 



(the softened gravitational n-body problem) are shad- 
owable for long times with as many as 150 phase-space 
dimensions (25 particles). We also demonstrate that, 
for the purpose of shadowing, a large n-body system 
can be viewed as a superposition of many small sys- 
tems, so that the distribution of shadow durations for 
a high-dimensional n-body system can be approximated 
using the distribution of shadow durations of many low- 
dimensional systems. 

The trajectories that we will attempt to shadow belong 
to a slightly simplified gravitational n-body problem in 
which there are N total particles, only M of which move. 
We do this because the shadowing algorithm we use takes 
time 0{M^) and we want to simulate a large system with 
a complex potential while keeping the time to compute 
a shadow tractable. Each moving particle interacts both 
with fixed particles and with other moving particles via 
Newton's gravitational force law, 
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where F^j is the force on particle i from particle j, 
and nij are their masses, rij is the distance between them, 
and £ is the gravitational softening parameter which, if 
non-zero, artificially smoothens the gravitational poten- 
tial in order to approximately emulate a system with 
more particles than are actually present and to avoid 
the singularity at = 0. We use normalized units 
in which each particle has mass the system has di- 
ameter of order unity, and the crossing time (the average 
time for a particle to cross the system from one side to the 
other) is of order unity We use a variable-order, variable 
timestep integrator for all integrations. We gener- 
ate noisy trajectories with local errors of about 10~^ per 
crossing time. To find shadows, we use an algorithm de- 
scribed in pO[ |, optimized to run between two and three 
orders of magnitude faster. Called iterative refinement, 
we use the same integrator as the noisy trajectory with 
tighter tolerance to estimate the full phase-space vector 
of local errors of the noisy trajectory, and then use a 
Newton-like correction to refine the trajectory until it 
has local errors as small as possible. For simple systems, 
the errors of the refined trajectory can be as small as 
the machine precision (10~^^), but the minimum local 
error acheivable with refinement increases as the number 
of dimensions increases, due to numerical errors in com- 
puting the Newton corrections. A trajectory produced by 
refinement is called a numerical shadow. The existence of 
a numerical shadow is expected to indicate the existence 
of an exact shadow of comparable duration [Tc| ] . 

A system with parameters N = 100, M = 1, e — was 
first shadowed by |0 , who found that the single particle 
can be shadowed for a few tens of crossing times, and 
that glitches (the point beyond which a shadow cannot 
be found) tend to occur more frequently during close ap- 
proaches with other particles. Combined with the fact 
that close approaches occur as a stochastic process [ p^ , 
we hypothesized that glitches also occur as a stochastic 
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FIG. 1: Histogram of shadow durations of 1000 unsoftened 
gravitational n-body systems. Each system has 99 fixed par- 
ticles and one moving particle. Noisy orbits have local error 
of about 10"'' per crossing time; numerical shadows were re- 
quired to have a maximum local error no bigger than 10"^*. 
The horizontal axis is in crossing times; the vertical axis is 
the measured probability density. The distribution fits an 
exponential curve with a mean glitch rate of about 0.07 per 
crossing time, indicating that the moving particle encounters 
glitches as a Poisson process in an unsoftend system. 



process. Fig. ^ shows a histogram of shadow durations for 
1000 AI = 1 systems. The distribution is well-fit by an 
exponential curve, suggesting that glitches are encoun- 
tered as a Poisson process in an unsoftened system. This 
is moderately surprising because it says that the number 
of crossing times that have occured in the past has little 
influence on the occurance of a glitch in the future, and 
thus there is in principle no natural upper limit on the 
duration of a shadow. 

The robustness of our shadowing algorithm with in- 
creasing problem dimension was tested in two ways. 
First, we contrived a slightly nonlinear Hamiltonian sys- 
tem designed to be easily shadowable {i.e., almost hyper- 
bolic), and successfully shadowed it for 20-50 Lyapunov 
times on average while the dimension was increased from 
2 to 180. Second, we searched for shadows of gravita- 
tional systems similar to the above consisting of 99 -I- M 
particles in which the M moving particles interact only 
with the 99 fixed particles. Such a system is equivalent 
to superimposing M uncoupled single-particle systems, 
and it will encounter glitches as a Poisson process with 
an aggregate rate M times that of a single-particle sys- 
tem. This results in average shadow durations that scale 
as 1/M. This was in fact observed in our simulations, 
giving us confidence that our shadowing algorithm works 
as well for large M as it does for small M. 

Fig. H demonstrates that the scaling is still 1/M even 
if the AI moving particles interact. This is moderately 
surprising because it suggests that although particles in- 
teract in motion, they do not interact in glitching. A 
possible explanation is that if 1 < M <C A^, then the 
coupling between the AI moving particles is weak on aver- 
age and we can still view the system as the superposition 
of M single-particle sy stems [pO[. Furthermore, we note 
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FIG. 2: Average shadow durations in crossing times of an 
unsoftened gravitational n-body system in which there are 
100 particles, Ai of which move, as a function of M. The 
noisy orbits have local error of 10~^ per crossing time; each 
numerical shadow was required to have a maximum error of 
lO"^'^. The points represent sample shadow durations, 30 
samples each for M = 1, 2, . . . , 19, 20, 25 and 10 samples each 
for M = 30, 35, . . . , 50. The "coupled average" line joins their 
averages, while 55/M°-^ and 55/Af^-^ are plotted for compar- 
ison. The "coupled average" line is statistically indistinguish- 
able from the "uncoupled average" one in which the gravita- 
tional interaction between moving particles is deleted, which 
both validates the robustness of our shadowing algorithm for 
large Ad, and suggests that even coupled particles encounter 
glitches independently of one another. 
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FIG. 3: Histogram of shadow durations of 250 systems iden- 
tical to those in Fig. ^ except with softening e = 0.1. Note 
that, in contrast to Fig. |^, the x axis extends to 1000 crossing 
times, and that the histogram height is zero near a shadow 
duration of zero, meaning that no particles undergo glitches 
until several crossing times have occured; in fact, of the 250 
systems sampled, the shortest shadow lasts 37 crossing times. 



that the cross-section for close approaches is not altered 
as M increases (which simply changes fixed particles into 
moving ones) , so the argument still holds independent of 
whether M < iV. 

With non-zero softening, [|l0| found that shadow du- 
rations increased significantly, while the correllation be- 
tween glitches and close approaches decreased. Fig. ^ 
shows the distribution of shadow durations for 250 M = 1 



FIG. 4: Average shadow durations in crossing times of a sys- 
tem with softening parameter e = 0.1 as a function of M. The 
points represent sample shadow durations, 9 samples each for 
Al = 1-5, 10, 20, 25, while the "coupled average" line joins 
their averages. Following the assumption that particles en- 
counter glitches independently of one another, the "overlay 
average" is artificially constructed for each AI — 1, . . . , 25 by 
superimposing M samples chosen at random from Fig. ^ and 
taking the minimum shadow duration of those samples. The 
resulting graph resembles the "coupled average" graph reason- 
ably well, except for the surprising effect of under-estimating 
shadow duration of the real system by about (20±10)%. 



systems in which softening has been set to e = 0.1, which 
is approximately half the average inter-particle separa- 
tion. The differences from Fig. |l| are quite marked. First, 
the distribution is peaked near 100 crossing times and 
has a long tail going out to hundreds of crossing times. 
Second, even though the local errors of the noisy and 
shadow orbits are identical to those from Fig. |^, the av- 
erage shadow length has increased by more than an order 
of magnitude from 14 to 218. Although this is roughly 
equivalent to the increase in the Lyapunov timescale [ p^ , 
the distribution is far from exponential. In fact, the most 
striking difference from Fig. is that the distribution has 
a vanishingly small density near zero shadow duration, 
in striking contrast to a Poisson process. In other words, 
virtually no particles undergo glitches until several tens 
of crossing times have occurred. If this remains true even 
for M > 1, then shadowing of softened gravitational sys- 
tems would be feasible even for large M, because the 
trajectories of all particles in the simulation would re- 
main valid for many crossing times. This question is 
addressed in Fig. ^, where we plot the average shadow 
duration for softened systems as a function of M, along 
with the shadow duration predicted by superimposing M 
single-particle systems. We see that although the dura- 
tion of shadows for coupled systems decreases as M in- 
creases, they decrease much more slowly than 1/M, and 
appear to be leveUing off at about 50 crossing times. This 
is consistent with superimposing single-particle trajecto- 
ries all of which have shadows that last several tens of 
crossing times, although it is surprising that the shadow 
durations are slightly longer than that predicted by su- 
perimposing single-particle trajectories. Now, if particles 
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encounter glitches largely independently of one another, 
and shadow durations for softened systems are long for 
most particles, then we can hypothesize that reasonable 
statistical results may be acquired from long simulations 
of large softened systems as long as only a small number 
of particles have undergone glitches, and the statistics 
taken depend on large numbers of particles. Thus, for 
example, the global spatial distribution of matter in a 
simulated galaxy may be correct, but the number of es- 
capers from a simulated globular cluster may be incorrect 
if the stars that escape happen also to be the stars that 
underwent glitches before escaping. 

The difference between the shadow durations for soft- 
ened vs. unsoftened systems undoubtedly is related to 
fluctuating Lyapunov exponents as discussed in [Q, |[ |j . 
Although we have not measured Lyapunov exponents 
explicitly, we measured a related quantity, namely the 
expansion and contraction factors across a timestep of 
the vectors that span the locally expanding and con- 
tracting spaces. In a uniformly hyperbolic system, these 
factors will always be greater and less than 1, respec- 
tively. An event of "non-hyperbolicity" can be observed 
by looking for areas along the trajectory where direc- 
tions which were previously expanding over long periods 
instead start to contract, or vice versa. If we plot the ex- 
pansion and contraction amounts along a trajectory, we 
find that "non-hyperbolic" events correllate well with the 
occurence glitches |2^. If we plot an M = 1 particle or- 
bit in three dimensions, we also observe that these events 
loosely correllate with times when the particle's orbital 
geometry changes in an obvious way. We postulate that 
the locally expanding and contracting directions of a par- 



ticle in the system are closely related to the geometry of 
the particle's orbit, so that changing the geometry of the 
orbit can cause these local vectors to become inconsistent 
as time progresses. In an unsoftened system, the geom- 
etry of a particle's orbit can be suddenly and violently 
changed by a close encounter. In softened systems, how- 
ever, there is no precise, short-duration "event" which 
triggers non-hyperbolicity; instead, the geometry of the 
orbit of a particle changes slowly, so that many crossing 
times occur before a glitch is likely. This helps to explain 
Figs. and |. 

In conclusion, we postulate that there is no feasible 
integration accuracy which will produce long shadows in 
unsoftened gravitational n-body simulations. This does 
not necessarily mean that such simulations should not 
be trusted, only that shadowing may be too stringent a 
measure of error. In contrast, we believe that there does 
exist a feasible integration accuracy for which softened 
systems are shadowable for many crossing times even for 
large n. Given the stark dependence of shadow duration 
on the softening parameter, however, we suspect that un- 
til a better understanding is acquired, the shadowability 
of physical simulations in general will need to be decided 
on a case-by-case basis. 

The author's software is available on request. 
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